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ABSTRACT 


A plane, heated jet of a constant property fluid exiting 
into a plane, straight channel of the same fluid is analyzed. 
The effect of buoyancy on the velocity and temperature fields 
is investigated. Two jet orientations are considered; a 
transverse jet, which enters the stream perpendicular to the 
main flow, and a longitudinal jet, which enters the stream 
parallel to the main flow. 

Solutions for the temperature and velocity fields for low 
Reynolds number flow are obtained using a finite difference 
scheme. Results are presented for three different values of 
the buoyant force, holding other variables constant. For an 
isothermal jet, results for the velocity field are presented 
at a higher Reynolds number. 
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NOMENCLATURE 


List of 

Symbols 

a 

Thermal diffusivity 

c 

Specific heat 

Ci 

Coefficients of interpolating polynomial 

D 

'N 

Channel depth 


Width of jet 

Dj* 

Non-dimensional width of jet 

/ 

F 

Unknown general function 

g 

Acceleration of gravity 

Gr 

Grashof number based on channel depth 

h 

Coefficient of heat transfer 

H 

Grid spacing in vertical direction 

J (A,B) 

Jacobian operator 

k 

Thermal conductivity 

Ki 

Coefficients of finite difference approximation 

M 

Square matrix 

P 

Pressure 

P (X/Y) 

Interpolating polynomial 

Pr 

Prandtl number 

q 

Heat flux 

Re 

Reynolds number based on channel depth 

T 

Temperature 

To 

Undisturbed stream temperature 

u 

Velocity in horizontal direction 
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u 

Non-dimensional velocity in horizontal direction 

u. 

D 

Non-dimensional jet velocity 


Average undisturbed stream velocity 


Surface velocity 

V 

Velocity in vertical direction 

V 

Non-dimensional velocity in vertical direction 

(x,y) 

Variable terms of interpolating polynomial 

X 

Space coordinate in horizontal direction 

X 

Non-dimensional space coordinate in horizontal 
direction 

y 

Space coordinate in vertical direction 

Y 

Non-dimensional space coordinate in vertical 
direction 

a 

Ratio of grid spacing in horizontal direction 
to spacing in vertical direction 

6 

Coefficient of thermal expansion 

e 

Non-dimensional temperature 


th 

Solution for 0 at i iteration 

M 

Absolute viscosity 

V 

Kinematic viscosity 

p 

Density 


Stream function 

Ip* 

Non-dimensional stream function 

ip*^ 

Solution for at iteration 

,2 

Harmonic operator 


Biharmonic operator 

< > 

Row vector 

{ } 

Column vector 

C ] 

Square matrix 


Order of 


5 



List of Subscripts 


a 

D 

i 

j 

jL 

ju 

L 

max 

s 

T 

u 


Ambient 

Downstream boundary 
Incoming stream 
Jet 

Lower edge of jet 
Upper edge of jet 
Longitudinal jet 
Maximum 
Surface 

Transverse jet 
Upstream boundary 
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I. 


INTRODUCTION 


Every coiranercial power generation plant releases some 
waste heat to the environment. Water from a nearby stream, 
lake, ocean is normally used as a reservoir for this 
rejected heat. As the demand for power increases, more and 
more heat will be discharged into these bodies of water. 

The advent of nuclear power plants, which have a lower ther¬ 
mal efficiency than conventional fossil-fueled plants will 
tend to increase the amount of waste heat per unit power 
output. With the emergence of these trends, concern about 
the effect of these thermal discharges has begun to grow. 

The utilization of a body of water for cooling purposes 
is normally carried out by diverting a certain mass flow 
from the body of vs^ater, passing it through a heat exchanger, 
and returning it at an increased temperature. Hence, the 
temperature of the body of water, at least in the local 
sense, is increased. 

From the biological viewpoint, there is no doubt that 
raising the temperature of water can be harmful to some forms 
of aquatic life. For example, for some fish a seasonal 
temperature rise of as little as 2°C can be fatal. Aside 
from such obvious damage, raising the temperature of water 
decreases its capability to retain dissolved oxygen, increases 
the metabolic rate of fish, and affects their reproductive 
habits and ability to resist disease. In addition, an 
increase in environmental temperature is beneficial to some 
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potentially harmful forms of plant life and bacteria [l]. 
Although the immediate effect of such changes may not be 
obvious, the normal equilibrium of the ecological system 
has been disturbed, and the possibility of long range 
damage does exist. 

In order to ensure that temperature disturbances are 
within safe limits, the precise nature of the temperature 
disturbance must be known. The quantitative value of the 
increase in temperature above the undisturbed temperature 
of the given body of water as a function of position and 
time would be of significant value in analyzing the effect 
of a thermal discharge. Several studies related to this 
problem have been carried out. 

The overall effect of a thermal discharge on a given 
body of water is one method of approaching the problem. A 
method for predicting the equilibrium temperature of a 
cooling lake was described by Sefchovich [ 2 ]. A study 
relating plant parameters to the average upstream tempera¬ 
ture in a stream or river was developed by Nahavandi and 
Campisi [3]. 

An experimental investigation of the temperature profiles 
in the vicinity of a cooling plant located on Yankee River 
was presented by Merriam [4], Approximate analytical models 
were compared with various experimental temperature profiles 
downstream of power plants situated on rivers by Polk, 
Benedict, and Lahey [S]. 

A possible model of a thermal discharge is a heated jet 
exiting into a given environment. For the case of thermal 
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discharge into a lake or ocean this model takes the form of 
a heated jet exiting into a quiescent fluid. An analytical 
solution for the two dimensional temperature and velocity 
profiles for the case of a heated, vertical , laminar jet 
exiting into an infinite environment was developed by Brand 
and Lahey [6], Additional work on a heated jet exiting 
into a quiescent fluid is now being conducted at Oregon 
State University [7]. 

The case of a heated discharge into a stream or river 
may be modeled by a heated jet exiting into a moving stream. 
An experimental investigation of a buoyant, turbulent, 
circular air jet exiting into a cross wind was conducted by 
Ramsey and Goldstein [8]. An analytical study of the 
surface spread of a submerged, buoyant incompressible jet 
was made by Tulin and Schwartz [9]. 

The problem undertaken in this study was that of a 
heated jet in a stream. This particular problem can be 
approached in several ways. The general problem is a tran¬ 
sient problem in which the temperature and velocity fields 
vary in three dimensions. Factors which must be considered 
are diurnal temperature and heat flux variations, tidal 
effects, seasonal temperature and flow variations, atmos¬ 
pheric conditions, stream meander path, bottom contour 
variations, fluid property variations, the possible turbulent 
nature of the effluent, and temperature stratification 
effects. Because of the complexity of this problem, several 
simplifications were made. 
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The transient conditions were assumed to be "slow" 


compared to the response time of the stream. Therefore, 
at any instant in time the stream could be considered to 
be at a quasi-equilibrium state. The time dependent terms 
were therefore not considered. The fluid properties were 
assumed to be constant, for simplification, and the axis of 
the river was assumed to be straight. 

It was decided to attempt to isolate the effects of 
buoyancy on the temperature and flow patterns. Also, to 
put the problem in a more tractable form, it was decided to 
consider variations in two dimensions only. Assuming that 
one of the dimensions would be the downstream dimension, 
there are two possible choices: a plane parallel to the 
river surface or a vertical plane passing through the axis 
of the river. Since the buoyancy forces act in the vertical 
plane, the vertical plane was chosen. In order to isolate 
buoyant effects, the bottom of the channel was assumed to be 
straight and parallel to the river surface. Also, for sim¬ 
plification, the river and the jet were assumed to be 
laminar in nature. The desired solutions were the tempera¬ 
ture and velocity profiles as a function of position and 
significant non-dimensional parameters. 
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II. ANALYSIS 


A. COORDINATE SYSTEM 

The problem under consideration in this work is flow 
in a channel into which a known quantity of fluid is injected. 
The injected fluid is at a higher temperature than the fluid 
in the channel. Since it is desired to study the effects of 
buoyancy, variations in the vertical plane only are consid¬ 
ered. The axis of the river is assumed to be straight. 
Rectangular cartesian coordinates are used, as shown in 
Figure 1. 

Two different orientations of the injected mass with 
respect to the main flow are considered. The longitudinal 
jet is injected parallel to the direction of main flow. 

The transverse jet is injected perpendicular to the direction 
of main flow. For the longitudinal jet, the origin of the 
coordinate system is located at the channel bottom, immedi¬ 
ately below the jet location. For the transverse jet, the 
origin of the coordinate system is located far enough 
upstream of the jet so that the incoming temperature and 
velocity profiles are not disturbed (see Figure 1). 

B. GOVERNING EQUATIONS 

Applying the Boussinesq approximation, which assumes 
that all density variations other than those giving rise to 
the buoyant forces may be neglected, and neglecting viscous 
dissipation, the steady two dimensional forms of the conti¬ 
nuity, Navier-Stokes, and energy equations are: 
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9x2 

9y2 


p[u + 

V ilj 

= - l£ + 


+ + pgSCT-T^] 

(3) 

9X 
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V 11] 
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(4) 

9x 

9y 

9x2 

9y2 



where T^ 

is the 

incoming 

stream 

temperature. 



By cross-differentiating with respect to x and y, and 
by subtracting equation (3) from equation (2), the pressure 
p is eliminated. The resulting equation is: 


u + V - U - V ^ 

8x9y Sx'^Sy 3y^ Sx'^ 9x9y^ 9x 


9x9y 9y^ 




(5) 


The continuity equation may be eliminated by defining 
a stream function , such that 


u = 

94) 

9y 

, V = 

94) 

9X 



Only two unknown quantities remain 

. to be 

determined, 

4' and T. The 

two 

remaining 

equations 

are 


ail* 

9y *- 9x 

94) 

9x 

9y 

= vV 4' ~ 

9T 

9^ 

(6) 

94) 9T _ 8iti 8T 
9y 9x 9x 9y 

= 

a V^T 



(7) 
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Equation (6) is the momentum equation expressed as the 
transport of vorticity, where the vorticity is defined as 

V . 

The equations will be non-dimensionalized by introducing 
the following variables: 


U = 


V = 


V. 


X 


X 

D 



where is the average undisturbed velocity, D is the 
channel depth, and Tj is the jet temperature. 

By substituting these into equations (6) and (7), 
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9V i{)* 


aifi* 

— [-^ 

ax aY 



Gr a 0 
Re2 ax 


( 8 ) 


le. _ Ml il = ■ 1 v2e 

aY ax ax aY 


RePr 
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where Re is the Reynolds number based on the channel depth, 

Pr is the Prandtl number of the fluid, and Gr is the Grashof 
number based on the channel depth and the temperature differ¬ 
ence between the jet temperature and the incoming stream 
temperature. 

C. BOUNDARY CONDITIONS 

Examination of the governing equations reveals that 
twelve boundary conditions are necessary to completely 
specify the problem. The necessary conditions are four 
conditions for i|) on x, four for i|) on y, two for T on x and 
two for T on y. These conditions will be specified on the 
channel bottom, the river surface, upstream and far down¬ 
stream of the jet. The boundary conditions for ip are first 
expressed in terms of velocity, and then related to the 

stream function from the definition of ip, i.e., = u, 

3y 

dip 

3x 

At the upstream boundary, the undisturbed incoming 
horizontal velocity profile is assumed to be a 20% truncated 
parabola, which is a common assumption in stream flow [lO]. 

A 20% truncated parabolic flow profile has the maximum 
velocity located two tenths of the total depth below the 
free surface (see Figure 2). For the longitudinal jet, 
this profile is modified by assuming that the stream 
velocity is equal to the jet velocity across the jet width 
(see Figure 3). For the transverse jet, the upstream 
boundary is located far enough upstream so that the velocity 
profile is equal to the undisturbed velocity profile. .The 
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Figure 2 - 20% truncated parabolic velocity profile. 
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Figure 3 - Velocity profile at upstream boundary 
for longitudinal jet. 
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vertical component of the stream velocity is assumed to be 
zero. For the longitudinal jet, the stream temperature is 
assumed to be equal to Tq except across the jet width, 
where the temperature is equal to the jet temperature, Tj. 
For the transverse jet, the upstream boundary is located 
far enough upstream so that the incoming temperature 
profile is also undisturbed. 

On the channel bottom, the no slip boundary condition 
requires that u = 0. The vertical component of velocity 
is also equal to zero, except at the transverse jet. In 
that case, v is equal to zero along the river bottom except 
across the jet width, where v is equal to the jet velocity. 
The river bottom is assumed to be adiabatic, yielding a 
temperature gradient boundary condition. 

On the channel' surface, it is assumed that there are 
no waves; i.e., v = 0. By assuming that the pressure is 
constant, and that the river surface is level, the appli¬ 
cation of the Bernoulli equation along the channel surface 
leads to the result that the surface velocity Ug is 

constant. Since U is constant, and the value of U_ at 

s “ 

the upstream boundary is known, the surface velocity is 
equal to the undisturbed surface velocity. At the surface, 
the conductive, convective, and evaporative nodes of heat 
transfer to the atmosphere are lumped into an overall heat 
transfer coefficient. It is assumed that radiative heat 
fluxes can be modeled by a known constant heat flux q. 
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A simple heat balance shows that the heat conducted into 
the surface plus the radiative heat flux into the surface 
equals the heat convected out. 

The velocity profile at the downstream boundary is 
assumed to be parabolic in form such that the total mass 
flow is equal to the undisturbed mass flow plus the mass 
flow of the jet. The vertical component of velocity is 
assumed to be equal to zero. It is also assumed that all 
heat added by the jet is convected out of the channel 
through the surface. Hence, the temperature profile 
returns to its original constant value T^. The downstream 
boundary is situated far enough downstream so that the 
above conditions are met within a given tolerance. 

The mathematical forms of these conditions are; 



U 


2 


x=0 


3 


£[8{Z)- 5(1) ] 



D D 


U 


j 


yjL<y<YjU 



yju<yiD 



V 






yjL<y<yju 
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yjU<y<D 


Trn = T 
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3y 
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Vt ~ 

- = 0 0<X<X-:l 

9x ■' 

= Uj 

II 

= 0 Xju<x 

^ 0 

ay 


9T^ 

W~ ' 

: 0 0lX<XjL 

Tt = 

= Tj 

9Tt 

= 0 Xj^<x 

3y 


y = D u = 

= iiiL = u„ 

9y 

V : 

= - ^ = 0 

Tx 

q • 

- k ^ = h(T-Ta) 

9y 

X->00 u • 

= il = Up (y) 

V 

= - ii = 0 

9X 

T 

= To 
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where the L subscript refers to the longitudinal jet, the 
T subscript refers to the transverse jet, the i subscript 
refers to the incoming velocity profile, the D subscript 
refers to the velocity profile at the downstream boundary, 
yjj^ is the y-location of the lower edge of the jet, yjy 
is the y-location of the upper edge of the jet, Xjj^ and 
Xjy are defined similarly for the x-location of the jet, 
and T is the ambient air temperature. 

Since the problem is cast in terms of the stream 
function ip, it is desired to reduce the boundary conditions 
to a more convenient form. At the undisturbed upstream 
boundary, the 20% truncated parabola gives the following 
form for the velocity profile. 

u = |i = [8 y - 5 (Z)^] 

ay 3 D D 

Integrating from 0 to D, it was found that the average 
velocity Uq is given by 



Now, integrating u (y) from zero to y, assuming 

i|;(x,0) = 0 

U^D « - 

rp = — [12 (Z)2 _ 5 (Z)3] 

^U 7 D D 

This is the upstream boundary condition for ^ for the 
transverse jet. For the longitudinal jet, if it is assumed 
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that the jet diameter is small with respect to the river 
depth, the expression for ii( becomes: 

'f' = 

= '/'u(y) + yjL<y<yju 

= 'J'u(y) Yju^yi'^ 

where Dj is the jet width, is the location of the lower 

edge of the jet, and yj^ is the location of the upper edge 
of the jet. 

For the longitudinal jet, the boundary condition at 
the channel bottom may be satisfied by assuming 41 is con¬ 
stant. Hence, 4 ) is chosen to be zero on the channel 
bottom for this case. 

For the transverse jet, the boundary condition for v 
may be integrated, giving 


4 ^.j. = 0 for 

X—X • 

4't = -UjDj (... JL' .) XjL<x<Xju 

4'T = "’^j^j ^>^ju 

The boundary condition for v at the surface may also 
be satisfied by choosing 4 ; = constant. The value of the 
constant is equal to the value of 4 ; at the surface as 
given by the upstream boundary condition. 

The downstream boundary condition for 4; may be found by 
integrating the downstream velocity profile, and evaluating 
constants from the surface and bottom values for 4 <. 
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After performing the indicated operations, and non- 
dimensionalizing, the final form of the boundary conditions 
becomes 


X=0 


i ( 12 Y^ - 5 Y-^) 


>U*OlY<YjL 

UjD.*(l^) YjL<Y<Yju 
Dj* 

^^,J;U*{Y) + UjDj* Yju<Y<l 


ax 
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Yju<Y<l 
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/ -u.D.*(- 2k) 
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3Y 


= 0 


OlX<XjL 


©T = 1 


XjL<X<Xju 


30. 


3Y 


= 0 


x.u<x 


Y=1 


3;|;* 

3Y 


£ 

7 


= 1 


'i’h* = 1 + UjDj* 


H - hD g ^ hD (Q _ j 
3Y K K A 


X-^~ il^L* = <j^D* (Y) = 


(12 + 3UjDj*)Y2 + (2UjDj* - 


'^T* = ’i'D*{Y) - UjDj' 


34/’* 

8X 


= 0 


0 = 0 


where the L and T subscripts have the same meaning as 
defined previously. 
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III. METHOD OF SOLUTION 


A. GOVERNING EQUATIONS 

The governing equations for the problem under considera¬ 
tion are: 




( 10 ) 


8i{)* d e _ 9i);* 96 

9Y 9X 9X 9Y 


1 


( 11 ) 


RePr 


In order to solve this problem using a finite differ¬ 
ence method, a finite difference operator for each 
differential operator in the governing equation is required. 
The method used to determine these operators is basically 
the same as that used by Dias [ll]. The operators developed 
in this study, however, are applicable to rectangular grid 
networks as well as to square grid networks. 

The procedure for determining the finite difference 
operators is to first define an interpolating polynomial 
P (x, y) having n undetermined coefficients. This poly¬ 
nomial can be used as an approximate to any general function 
F. (F can be interpreted as a stream function, for example.) 
By requiring that P (x, y) be equal to F at each of n 
points, n equations involving the n undetermined coeffi¬ 
cients are generated. Hence, the coefficients of P (x, y) 
can be determined in terms of the value of F at each node. 
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By applying the desired differential operator to P (x, y), 
an approximation for the operator in terms of the value of 
F at each node is obtained. 

To develop the biharmonic operator, a polynomial must 
be chosen such that the fourth partial derivatives with 
respect to both x and y are non-zero. The grid network 
through which P (x, y) is to be passed is shown in Figure 
4. Since there are 25 points in this grid, 25 unknown 
coefficients are required in the polynomial P (x, y). 

Hence, an appropriate choice for P (x, y) is 

P(x,y) = C^y^x'^ + C 2 y^x^ + C 3 y^x^ + C 4 yx^ + C 5 X^ 

+ Cgy'^x^ + C^y^x^ + Cgy^x^ + C^yx^ +0^0^^ 
tC^iY^x^ +^247x2 +C^^x^ 

+Ci6y^x +Ci7y3x tCigy^x tC^^gyx +C 20 X 

+C 2 iy'^ +C 22 y^ ■'■^23^^ +C24y +C 25 ( 12 ) 

This equation may be written in matrix form as: 

P{x,y) = <Vi(x,y)> {C^} i = 1,2,...,25 (13) 

where (x, y) is x^y*^ corresponding to in equation (12) . 

This polynomial is now required to pass through the 
value of Fj at each node j, where Fj is the value of the 
function F at node j. At each node j, the following 
relation is obtained: 

P(Xj,yj) = <Vi(Xj,yj)> {C^} = Fj i = 1,2,...,25 

(14) 
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By applying this relation at each of the 25 grid points, 


the following set of equations results; 


\ 

1 

A 

>1 

r—1 

X 

> 

V 

1_ 



^0 1 

r 

<Vi(x2 y2)> 

• 

< 

C2 

• \ 

V ^25/ 

• 

<Vi- (x^ ,y ) > 

1 25'^25 

U- »J 

1 



(15) 

j = 1,2,...,25 


Or, in matrix notation: 


{Fj} = [M]{Cj} 


j = 1,2,...,25 


(16) 


Solving this matrix equation for Cj, 


{Cj} = [M]~^{Fj} 


j = 1,2,... ,25 


The coefficients Cj of P (x, y) are now known, and 
P (x, y) can be used as an approximation for F. The 
biharmonic operator may now be applied to the polynomial 
P (x, y). 

v4(P(x,y)) = <v4 ( v^U,y)> {C^}) 

= v4<(Vi(x,y))> (Ci) i = 1,2,...,25 

Since an operator for the biharmonic operator for the 
central node is desired, this expression is evaluated at 

xi3,yi3. 

v4(p(x,y)) = <v4 (Vi(x,y))> [m]~ 1 {F.} i=l,2,...,25 

^ (17) 

= <Ki> (Fi) i=l,2,...,25 (18) 
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where <Kj^> is a row vector of constants obtained by 
multiplying the matrix expressions appearing in equation 
(17). The result is that are the coefficients of a 
finite difference approximation for the biharmonic operator 
applied at the central node of the grid in Figure 4. 

In order to ensure geometric consistency, this same 
interpolating polynomial was used to develop the operators 
appearing on the left hand side of equation (10). The 
desired operator is applied to the polynomial, and evalu¬ 
ated at the central node. The results for a = 5 are 
presented in Figures 5, 6, and 7, where a is the ratio of 
the spacing in the x-direction to the spacing in the 
y-direction. 

A finite difference approximation for the harmonic 
operator, which is required in equation (11) is available 
in the literature [12], Consistent approximation for the 
other operators appearing in this equation are also 
available. These operators are presented in Figures 8 
and 9. 

B. BOUNDARY CONDITIONS IN FINITE DIFFERENCE FORM 

Assuming there is no jet present, the boundary condi¬ 
tions for ip* are all of the form 

il>* = f(t) 



where t is the space coordinate tangent to the boundary. 
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irence approximation for V operator for uniform 
grid (a=5), applicable at central node. 
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rectangular grid (a=5), applicable at central node. 
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Common Factor 



Truncation Error = 


Figure 8 - Finite difference approximation for 

operator applicable at central node 



Common Factor 


2H 


Truncation Error = 


2 


) 


Figure 9 - Finite difference approximation for 

9/3X or 3/3Y applicable at central node 
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n is the coordinate normal to the boundary, and f (t) and 
g (t) are the values specified at each particular boundary. 
The first condition on \p* is satisfied in the finite 
difference scheme by evaluating f (t) at each node on the 
given boundary. These values are then used in the finite 
difference computations. 

To satisfy the boundary condition for the normal deriva¬ 
tive of rp, g (t) is evaluated at each node on the boundary. 

A finite difference approximation for the normal derivative 
is then required. An offset three point operator is used 
because the truncation error is consistent with the trunca¬ 
tion error of the governing equation. The finite difference 
operators used at the boundaries are presented in Figures 
10 and 11. 

The temperature boundary conditions assume various 
forms. At the upstream and downstream boundaries, the value 
of the temperature at the boundary is specified. At nodes 
located on these boundaries, the known value of the temper¬ 
ature is used in the finite difference computations. On 
the channel bottom, the derivative of the temperature is 
specified. At nodes on this boundary, the offset three 
point operator utilized previously is applied. At the 
channel surface, the sura of the temperature and its deriva¬ 
tive is specified. At nodes on the channel surface, a 
mixed boundary condition exists. In this case, the value 
of the temperature at the node located on the surface is 
added to the three point operator shown in Figure 11. 
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For most computations, the jet center is located between 
two nodes. In addition, the jet diameter is assumed to be 
less than the grid spacing in all cases. The effective 
diameter of the jet is equal to the grid spacing on the 
boundary on which the jet is located. This is a result of 
the fact that any length smaller than the grid spacing 
cannot be accurately modeled in a finite difference network. 
This can best be illustrated by example. Figure 12 depicts 
two possible jet configurations, both having the same volume 
flow rate. Also shown are the stream function profiles, 
superimposed on a grid network. In both cases, the jet 
width is less than the grid spacing. It can be readily 
seen that the values of the stream function at the two 
nodes on either side of the jet are the same, regardless 
of the jet diameter. The value of the stream functin is 
changed by the amount of the volume flow rate of the jet. 

If the jet width is smaller than the grid spacing, the 
effect of the jet on the finite difference grid is a fixed 
flow rate between two grid points. The effective diameter 
of the jet is then equal to the grid spacing, and the 
effective jet velocity is the flovr rate divided by the grid 
spacing. 

By modeling the jet in this manner, the edges of the 
jet fall on nodes of the finite difference grid. In a 
continuous system, there is a step change in both the 
temperature profile and the velocity profile at these 
points. This step change cannot be modeled precisely in a 
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Figure 12 - Velocity profiles and resulting stream 
function profiles for two jets located 
between nodes 1 and 2. 
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finite difference scheme. By specifying the temperature 
at these two nodes in the finite difference grid, the 
temperature profile is as shown in Figure 13. It can be 
seen that the total area under the temperature profile 
curve is equal to the area under a step change of tempera¬ 
ture 2 X Tg. The energy input of the profile shown in 
Figure 13 is equal to that of the step change of height 
2 X Tg. This temperature may be defined as the effective 
jet temperature. 

If the three point finite difference operator shown in 

Figure 9 is used to calculate the velocity profile based 

on the relation V = , it can be shown that the effec- 

3Y 

tive jet velocity defined earlier is consistent with the 
definition of the effective jet temperature. 

The above examples are directly applicable to the 
transverse jet. By analogy, the results may easily be 
extended to the longitudinal jet. 

For the results used to construct the convergence plot, 
the jet center is located at a node. For this orientation, 
the effective jet diameter is twice the grid spacing. 

C. NUMERICAL SOLUTION SCHEME 

In the numerical solution of this problem, the governing 
equations are rearranged 

= - 2£ 11 + [111 3 (18) 

Re 3X 3Y 9X .9X 9Y 

= RePr [11_ 11 - H] (19) 

9Y 9X 9X 9Y 
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Figure 13 - Temperature profile over finite difference 
grid in the vicinity of the transverse jet. 


39 
















These equations may be expressed in an alternate form if 
we recognize that the Jacobian operator is defined as 


J(A,B) 


8A 8B 3A 8B 
dX 8Y 9Y 9X 


Equations (1) and (2) may then be rewritten as: 


= - H H - J('i'*,V^i|;*)Re (20) 

7^0 = RePr J(ti;*,0) (21) 

In the solution of these equations, an iterative scheme 
based on the uncoupled linear homogeneous forms of these 
equations is used. Once this linear system of equations is 
solved, the right hand side of equations (20) and (21) may be 
calculated. 

The linear equations are then resolved, using the calcu¬ 
lated values of the right hand side of the equations as known 
inhomogeneities. This procedure is repeated until the desired 
convergence criteria is satisfied. 

The uncoupled linear homogeneous forms of equations (20) 
and (21) are: 

= 0 
7^0 = 0 

The method of solving these equations with the boundary 
conditions for the transverse will be illustrated using the 
grid in Figure 14 as an example. 
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The boundary conditions on specify the value of i|)* on 
all boundaries. Therefore, the value of tj;* at nodes located 
on the boundaries is known. The value of ip* at each interior 
node is unknown. In order to solve for the value of ip* at 
interior nodes, one equation for each node is required. 

For all ip* values not adjacent to any boundary (nodes 15, 
16, 21, and 22 in example), an equation is obtained by applying 
the governing equation at the node. The finite difference 
operator shown in Figure 5 is used. 

For all \p* values located at nodes adjacent to the upstream 
boundary (nodes 8-11 in example), an equation is obtained by 
applying the boundary condition 14— =0 at the adjacent 
boundary node (nodes 2-5 in example). The finite difference 
operator shown in Figure 10 is used. 

For <p* values at nodes adjacent to the downstream boundary 
(nodes 26-29 in example), an equation is obtained by applying 
the boundary condition 14" “0 the adjacent boundary 

node (nodes 32-35 in example). The finite difference operator 
shown in Figure 11 is used. 

For ij;* values at nodes adjacent to the channel bottom, but 
not adjacent to either the upstream or downstream boundary 
(nodes 12, 20 in example), an equation is obtained by applying 
the boundary condition =0 at the adjacent boundary node. 

The finite difference operator shown in Figure 10 is used. 

For \p* values at nodes adjacent to the channel surface, but 
not adjacent to either the upstream or downstream boundary, 
(nodes 17, 23 in example), an equation is obtained by applying 
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9 it* 9 

the boundary condition = sj- at the adjacent boundary 

node (nodes 13, 19 in example). The finite difference operator 
shown in Figure 11 is used. 

Since the value of it* is known at nodes located on the 
boundary, whenever a finite difference operator involves the 
value of it* on a boundary, the known value of it* is used. 

The resulting system of equations is a set of algebraic equa¬ 
tions for each unknown. This system of equations is solved 
by Gaussian elimination. 

The boundary conditions for 0 specify the value of 0 at 
the upstream and downstream boundaries, and at the two nodes 
on either side of the jet. However, the value of the tempera¬ 
ture on the channel surface is unknown, and the value of the 
temperature on the channel bottom, except for the nodes on 
either side of the jet, is also unknown. Again, one equation 
for each unknown is required. 

For each 0 unknown not located on a boundary (nodes 8-11, 
14-17, 20-23, 26-29 in example), the equation V^0 =0 is 
applied at the node at which the unknown is located. The 
finite difference operator shown in Figure 8 is used. 

For 0 unknowns located at nodes on the channel bottom 
(nodes 19-25 in example), an equation is obtained by applying 
the boundary condition =0 at the node. The finite 

difference operator shown in Figure 10 is used. 

For 0 unknowns located at the channel surface, (nodes 
12, 18, 24, 30 in example), an equation is obtained by 
applying the boundary conditions ^+^0=0 at the node. 
The finite difference operator shown in Figure 11 is used-. 
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The result is an equation for each 6 unknown. This system 
of equations is solved using Gaussian elimination methods. 

The method of solving equations (20) and (21) is presented 

in the form of a flow chart in Figure 15. The superscripts, 

£ th 

e.g., ip* refer to the solutions at every node on the Z 

iteration in a particular loop. The quantity Rg (i/)*"', 6^) is 

defined as the right hand side of equation (20) based on the 

known values of and 0^. The quantity Rg is the right hand 

side of equation (21) defined similarly. 

The convergence tests for ip* and 9 are: 


,i-l ,i 




< 0.001 






< 0.001 


It is assumed that when the above criteria are satisfied 
at every node, the solution has converged. 

It should be pointed out that the governing equation is 
not applied at every unknown node. The non-linear terms 
calculated from the values of functions from the previous 
iteration affect the known value in the governing equation only 
at nodes at which the governing equation is applied. However, 
the solution for every node is affected, because the equations 
are solved simultaneously by Gaussian elimination techniques. 

In the calculation of the right hand sides of equations (20) 
and (21), two different techniques were used. The first 
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SOLUTION 


Figure 15 - Flow chart for iterative scheme, 
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technique utilized the operators shown in Figure 6, 7, and 9, 
which were derived earlier. 

A second method , which provided improved numerical 
stability, was also used. This method is described in [13]-. 
This is simply an alternate way of calculating the Jacobian, 
which appears both in the right hand side of equation (20) 
and equation (21). 

Using this method, the Jacobian Operator is expressed in 
three alternate forms: 


J(A,B) 


dA 9B aA 3B 
3X 3Y 3Y 3X 


3 

3X 

3 

3Y 


(A 

3B. 

3y^ 

3 

3Y 

(A 

3B. 

3X^ 

(B 

3A. 

3X' 

3 

3X 

(B 

3A. 

3Y^ 


When the quantities are grouped in this manner, applying 

the three point central difference operator shown in Figure 9 

leads to three different numerical methods of calculating 

J(A,B). By calculating J(A,B) by each of these methods, and 

using the average of the three values, some forms of numerical 

instability are eliminated. 

This method may be applied directly to equation (20). 

However, in order to apply it to equation (21), the value of 
2 

7 tp* must be calculated at each node. This is done by applying 

the finite difference operator shown in Figure 8. The method 

2 

described above can then be applied to calculate J(ij<*,V ifi*) . 

The computer program used to obatin the results which are 
presented is included in Appendix B. 
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IV. RESULTS 


A. PRESENTATION OF RESULTS 

In the results, velocity data are presented in the form 
of streamlines in the channel. Temperature data are presented 
in the form of isotherms in the channel. Both the streamlines 
and the isotherms are presented in non-dimensional form, and 
are plotted with respect to non-dimensional values X and Y. 
These graphs are based on linear interpolation between calcu¬ 
lated values of i/)* and 0 at nodes in the grid. It should be 
noted that in the graphical plots for both and 0, the scale 
in the Y-direction is twice that in the X-direction. The grid 
used for all calculations consisted of 7 nodes in the Y- 
direction and 19 in the X-direction. 

For all results presented, the surface boundary condition 
on the temperature is ^ ^ 0 = 0 . Additionally, the 

ratio of the flow rate of the jet to the flow rate of the 
river is 0.2 in all cases. This value was chosen as a repre¬ 
sentative value for an actual power plant, and was felt to be 
small enough so that the surface boundary condition on the 
vertical velocity (i.e., no waves) remains valid. 

In Figures 16-21, the streamlines and isotherms for a 
transverse jet of width 5D/6 are presented. The Reynolds 
number and Prandtl for all these figures are Re = 0.3 and 
Pr = 6.0. Results for three different values of the ratio 
of the Grashof number to the Reynolds number are presented; 
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Gr/Re = 600, Gr/Re = 0, and Gr/Re = -300. The case of a 
negative value of Gr/Re implies that the temperature of the 
jet is less than the temperature of the stream. Hence, the 
buoyant force acts in the opposite direction, i.e., downward. 

In Figures 22-27, the streamlines and isotherms for a 
longitudinal jet of width of D/6 located between Y = 1/3 
and Y = 1/2 are presented. For all these figures. Re = 0.3 
and Pr = 6.0. The three values of Gr/Re presented are 
Gr/Re = 600, Gr/Re = 0, and Gr/Re = -300. 

In Figures 28-29, the streamlines for both a transverse 
jet of width 5D/6 and a longitudinal jet of width D/6 
between Y = 1/3 and Y = 1/2 are presented. The values of 
the non-dimensional parameters for these figures Re = 40, 

Pr = 6.0, and Gr/Re =0. No isotherms are presented because 
the jet is assumed to be at the same temperature as the 
river. 

In Figure 30 and Table 1, the surface temperature as a 
function of X is presented. The temperatures for the longi¬ 
tudinal jet are presented in tabular form becasue differences 
in the temperatures for the three different values of Gr/Re 
are difficult to distinguish graphically. For all cases in 
Figure 30 and Table 1 the Reynolds number and the Prandtl 
number are Re = 0.3 and Pr = 6.0. The value of Gr/Re for 
each case is indicated in the results. 

B. DISCUSSION OF RESULTS 

The effect of the buoyant force on the flow patterns and 
isotherms can be seen by examining the graphical results in 
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Figure 30 - Surface temperature vs. X for a transverse jet 
of width 5D/6, Re = 0.3, Pr = 6.0. 
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TABLE I 


Surface Temperature {0_) 

W 

Longitudinal Jet 


X 

Gr/Re = 600 

6 

s 

Gr/Re = 0 

Gr/Re = 

0.0 

o 

• 

o 

0.0 

o 

• 

o 

0.833 

0.1083 

0.1077 

0.1074 

11.667 

0.0907 

0.0901 

0.0898 

2'. 500 

0.0632 

0.0635 

0.0636 

3.333 

0.0431 

0.0437 

0.0440 

4.167 

0.0294 

0.0300 

0.0304 

5.000 

0.0200 

0.0300 

0.0304 

5:. 833 

0.0137 

0.0141 

0.0144 

6.667 

0.0094 

0.0097 

0.0099 

7.500 

0.0064 

0.0066 

0.0067 

8.333 

0.0044 

0.0046 

0.0047 

9.167 

0.0030 

0.0031 

0.0032 

10.000 

0.0021 

0.0022 

0.0022 


-300 
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V 


•J. 


FigTon r e s 1®—ZL ZZ*—2277.. case- off G£’/Ree==OC implies the 


Imoyaiifc effect its. rre^iiij^iijfLhB iir. ttie: momenttimreqaation;.. In 
this case^ ncc bunYantt fbccam acts? on. thee flow/, andctheumomentum 
egraiaticiii is irndepeaideirtt off thees tTKnperatureeeqnatzonr.. A posi- 
ti\re value of Qr/Rte iTnp-Tvfyiy thcat:: the' tempesatureeof ithe jet is 
greater than T^.. liT. th±s case^ the buoyant-: force iaets in the 
pasi-txve I--<ifrect±'car.. Ai. nega2tive3 vaiueioff G£/Re=implies that 
the temperatuxs a£ thee jeCfct Uss less' tharr.T"^ . IQr.this =case, 
the baoyaiit foxce. actsv ixr thei negative Yrdirectionv.. In order 
to examine the effect: off Buoyaxcy , the:-casessihrwhich Gr/Re 
is non—xero axe campaxed'. toe the; case ie whichhGr/Re;=-0. 

For the txarrsverse jo-tv, the 'P* — 0 streamline;indicates 
the uppex boundary.- of- the flow introduced:by the: jet. No 
flow e-ver errosses thie boundary.. Because;theeflow. introduced 
by the jet occupies the space: below thisrboundary^ the flow 
originally in. the chaainel: ife flowing through:.a:.smaller area, 
and is therefore accelerated: downstream-off the: jetz . 

A positive value of'Gr/Re-causes the; transverse jet to 
penetrate further into the stream than the jet-:fdr:which 
Gr/Re = Q, particularly near the jet location.. AS=a result 
of this larger "obstruction" to the flow in the channel, local 
regions of decelerated and accelerated flow are created near 
the jet (Fig. 17).. A. negative value of' Gr/Re - restricts the 
penetration, of the jet.. Dr this case, theeflow. originally 
in the channel is:- abl& toe surpass the. jet:in.-a:smooth manner 
(Fig. IS);. 
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gear tdtBe rocrr^iAaidinail tfaeeeffeet'-of the buoyant 

feuroa CEir. thee ffljocw'isr not" ass pronounced Jas .in the transverse 
Th±ff isE becausB tiie' initial 1 velocity of the jet is in 
ths- dawnrs trea nr diireclrion, andd therefore >. the heat added by 
t±rg ji'et: is ennveertexi' downstreanr mores quickly than for the 
tirangverses j^edr.. Hence , thes local-values. of 0 are not as 
great in. thes casse of: the' Ibngitudinall jet. Since the 
Emoyant force: is proportionail toe 0 > ,thesbuoyant forces near 
the j-et aECEs noct as great in" the-;-cases of: a longitudinal jet. 

The edEferdr of the buoyant-;force:can; be seen by comparing 
the stxeamrinee for the case: in: which Gr/Re is greater than 
zero to the case in which Gr/Res=-0;. The positive buoyant 
force causes" the streamlines-to:be'bent upward with respect 
to the case in. which Gr/Re —0' (Fig; .22}, This is due to 
the velocity in. the positive'Yrdirection caused by this 
b^layarrt farce.. A', negative value--of: Gf/Re has the opposite 
effect (E'ig^.. 233),, causing the: streamlines to be bent 
down-ward.. 

The effect, of' buoyancy on the temperature field for 
both jets can be seen by examining the isotherms presented 
in Figures 19—21 and 25-27. At downstream locations near 
the jet, temperatures are greater.for a positive buoyant 
force than for the case in which there is no buoyant force 
present.. This is-due to the; fact that-.a positive buoyant 
force con-tributes: to the heat; transfer; in the vertical 
directiorr by.- increasing the:vertical 1 velocity component. 
Hence,, the heat added by the jet is dispersed more quickly 
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im tixe ver±±c3Ell dxrectlort,. andi temperatures near the jet 
sirs higher:.. Tlriss results^ inr.anr.increased surface temperature. 
Bi£ca.usg aJjX. heMr added' by thee jet~inust~ leave the channel 
tthrrrugh. thes euxrfarre'^. and'thee amDunt"of heat transferred at 
the; surface: iis propox.tionall toe thee-surf ace temperature, more 
heat rs trarrsferred. out off thee channel near the jet for a 
posEiitiv.e'- buc3Yant for.ee than.eforraezero buoyant force. 
Hherefore thee temperatureseatedownstream locations are lower 
ferr the. poartive: buoyant force;. This effect can be seen by 
^amxaing thee surface temperature:data presented in Figure 
20.. This effect occurs for both:.the longitudinal and trans¬ 
verse: j-ets,, but is'more pronounced for the transverse jet. 

R. negative-, buoyant force has- the'opposite effect, causing 
the temperatures near the jet:to:be.smaller, and the down¬ 
stream temperatures higher. 
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W.. DISCUSSIDN' 


A. OI}C!i:7EE(^CE CHF ETNEAR- SOLUTION;'FDR? STREAM FUNCTION 

Tfoe valxdxty of thee niatxxx off linEarrequations for the 
stream farurtion. can be:- verified by constructing a conver¬ 
gence plat to: determines the value' to: which" the numerical 
EOlutiau WQU'ld; canv.ergie if^ an infinitee number of nodes were 
used. This value, may then be: compared: to:an. exact solution. 

A convergenee platt re constructed; byyplotting the value 
of the numeirical sdutron obtained at. a: given point in the 
field of the problem versus 1/N , where: N is the number of 
divisions per side of: the field. By making the number of 
divisions greater,, different values, of: the: numerical solution 
cire obtained at the same point in the-field. The plot of 
the value of the aalution. versus 1/N.^'may then be extrapo¬ 
lated to determine the point, at which: the-.plot .crosses the 
1 /n 2 = 0 axis.. This vafLue of the solution is then taken as 
the value to which: the' numerical solution:.would converge if 
an infinite number of divisions were possible. 

In order to apply this method to the linear (zeroth 
order) solution, for the stream function, a'fixed field must 
be chosen. The field chosen was a section of the channel 
of length 5D.. The case considered was: a, longitudinal jet 
of width 2H located at the upstream boundary and in the 
center of the channel... By utilizing: a: grid in which the 
spacing in the X'-direction is five'times--the spacing in. the 
Y-direction, a grid with an equal number of divisions on 
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eacfit edge of thie ffiifiCd:. possibles . CSnvergence plots were 
camstrxtcrted: at: the£ points: X = 2:. 5^ , Y;’== 0 ^ 25 and X = 2.5, 

Y = Q«T5* Tbna mesit: sl..zes~ thati wereeusedc'were N = 4, N = 8, 
and K = 12^ Tli& N' = 44 gxid' was- tooc coarse to yield useful 
convergence: data.. SHtce: a grid:, corresponding to N = 16 
requires excessive- computer, s.toragev , only the solutions 
correspondi’rrg to: n: == 8: and" N'— 121 werefused in the extrapo¬ 
lation process.. The convergence-pibts£ for these two points 
are shown in. Eigurers; 3X and 32.'., 

An exact sniut±an at these: two:, points is obtained by 
noting that the surface and bottonr boundary conditions for 
the velocity- are. the same as those: for' Couette flow between 
flat plates., The specification ofa: knov/n flow rate for 
incompressible flow is- equivalent tor specifying a driving 
pressure gradient in the downstream direction. Hence, for 
flow in -the channel far enough downstream so that the effect 
of the jet on -the flow has been damped, out, the exact 
solution is the solution for Couette'flow between flat 
plates with an imposed pressure gradient. This solution is 
given in Schlicting [14]. It should be noted, that this 
velocity profile is parabolic in y, and that U (0) = 0 and 
U (D) = Ug. For an incompressible flow, the pressure 

gradient IZ. may be determined from a. specified flow rate. 

8X 

At -the downs-tream boundary, several, assumptions are made. 
It is assumed' tha-t the: effect, of: the: jet .may be neglected, , 
wi-th the ®£cept±arr of the increased: flow added by the jet, 
and that v = 0.. It is also assumed that the velocity profile 
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Eignxe 31 — Convergence plot for stream function at 
X 2.5, y = 0.75. 
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Figure Z2 — Convergence plot for" stream function 
ai X = 0.25, Y = 2.5. 
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is pff.TnThfT.T:r(r. iirr; ^v, aoidi thee fibw.'. rateeisiknown. In 

ai2a£t±caT.„ thee ved!brd±x-'atr: y/= O' iss equal 1 toe zero and that 
ttfrffi v^DcnriitY" ah. y’ = EE’ iss eguall toclTg, . Sinceethe velocity 
BorcrgtT'g- ids pcErsEhallite iir. y.^, and' sincee thee three conditions 
which, d'eterraines thee CDefidjclents-of f the; parabola are 
iden.tical- ta tliee cnndritions' for Couette-flow with an imposed 
pressHur^ gr^tdient,, the velbcity profile-atethe downstream 
heur n dary is iderrtica'll tor the exactr solution; 

it was atssumed:*. that; the' exact-soiution.“was valid at 
S = 2..5.. The vaiiieg obtained from.-theeconvergence plot may 
then, be compared tor thoB.e:-obtained from: the exact solution, 
at the point- X'. = 2L.5c,. Y: == 0.25, the:-exact , solution is 
i|;* ==• Q..12r723:;' the sorlu-tion obtained, from- the convergence 
piot is 4** — 0[.JI2r73i.. A'tr the point" X'== 2 d 5 > Y = 0.75, the 
exact solution is- ii;*" = 0^83170; the. solution obtained from 
the convergence plixtr isr i(;*’ — 0.83163; . Asrarresult of the 
agreement-between, these vaiues, it"_may be • concluded that 
the linear solution for 4 ;*', which is the zeroth order solu¬ 
tion for the iterative scheme, is valid. 

B. NUMERICAL STABILITY OF ITERATIVE SCHEME 

By formulating the problem in terms of the linear 
solutions as shown in equations (20) and (21) , the nonlinear 
terms in the energy equation are multiplied by the Peclet 
number Ee,, where P.e =-Re-. Pr_. The; nonlinear: terms in the 
momentum equaticar. are multiplied by Re-; . However, using the 
iterative' procedure described earlier; theenumerical solution 
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£mr tdie eareEC try - eEcpxatlon beeomess unstable at approximately 
Fe = 2.. ffiiinnoes tliec energy• arrdimontentum'equations are coupled, 
frTT<» m T t-iirr r- sal]ut±an. becomess uirstableeat-Pe = 2. 

SfeweEcadl raethads^ were emplbyed:' inr.an- effort to eliminate 
tEdis xnrstniiiiriity.-.. The methiodi initially used to calculate 
tbe nonlinear' terras^ was to: use', thee operators shown in 
Figuree 5-,. 6l„ anx£ 7.. Usingr thisemethod, the energy equation 
feecame unrsEtable SEt approximately • Pee =-1'. When uncoupled 
fr o m the earergy equation, thee momentum equation became 
unstable ach Re = 15. 

The method' of calculating- the-; nonlinear terms in terms 
of the GTacobtarr. aES described: earlier; was then attempted. 

Using this method, the Pecletrnumber, at which instability 
oxDcurred in. the energy equation was.=^ raised from 1 to 2. 

The Reynolds number at which instability occurred in the 
momentum' equation ,, uncoupled: from-the: energy equation, was 
raised from 15 to 50. 

Oecreasinq the grid spacing*was-also attempted. This 
had no noticeable effect on the stability of the system of 
equations.. 

An attempt was made to begin with the stable solution 
for Pe = 2.0, and increase the Peclet number in steps of 
1/25. Using this procedure, instability occurred at Pe = 2.15. 

C» CONCLUSIDNS 

Eecau'se' all. solutions pr'esented; in this study are based 
on itexationg from the linear:^, uncoupled solutions of the 
derived governing equations, the validity of the numerical 
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usEsdi toe abtaiir. tiiesee solutions is of utmost impor- 


ttSHTCE:.. Tlte v-aiMiiity of:: tfaeeoperators eand method of solution 
fianr theB Hrrmair energy equationr;has sheen demonstrated by its 
freguent: xisee with: successs inr. theeliterature [12, 15], The 
itKsn'Xts cdC tiies crmvergence-'pibtr indicate that the linear 
ggluttan. of the: Linear momentum' equation is also valid. 

EsESEed: car theses ohrservations^ , it‘: is s felt that the solutions 
presented- iim thiss reportr areecalculated from an essentially 
SGrund basis.. 

The esEfect o±'buoyancy on:: thes flow patterns and tempera¬ 
ture distritu-tions has beenr demonstrated for the cases 
cenTsidered',. However, the: original 1 motivation for this 
effort was -toe demonstrate the: effect .of buoyancy on the flow 
in-, an. actuad physical channel!.. Forrsuch channels, the 
actual values of'Reynolds and'Peclet-numbers experienced 
ace; grea-ter;-tharr 10^. Because:ofithe:restraints on the 
range of; Rd and! P.e. due to nonlinear: numerical instability, 
it was not possible to generate:solutions which are appli¬ 
cable to actual channels. 

The maxiraxim value of the Peclet number for which results 
were obtained was approximately 2.0. Using the value 
Pr = 6 as a typical value of the Prandtl number for water, 
■the maximum Reynolds number for which results for the 
caup:Ied temperature and flow fields'were obtained was 
approx imately • 0!. 31.. 

fis a. result of these limitations> some of the assump¬ 
tions made for the boundary conditions for the longitudinal 
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jet nrs-y fee (^ueetrroned:.. Ijt, the fbmralationr.ofithe boundary 
condition: fenr thre temgemture' foxr thee longitudinal jet, it 
was assumed that ntx hesh wa^ cDnductediupstream“to the 
X—location of the jet.. Th±s assumption”, issvalid: for 
Pe >> L. Since the maximuir Peclet. nuinberrused: in the cou¬ 
pled profeLeir. wae Ee = this assumption” is5questionable. 

For the velocity bound'any condition:, forrthe : longitudinal 
jet, it was- assumed. thaEt: the distrubancee toe the:-flow field 
caused fey the jet does net: propagate-upstream.” to, the 
X-location of the j;et.. Hris assumption.”, issvalid for 
Re >> L. Because the maximum Reynolds:numbereused in the 
solution of. the coupied energy and momentum-equations was 
0.3, this assumptian is no longer valid. 

D. RECOMMENDATION S 

Using- a finite difference method: off solution-, one 
algebraic equaticor. is generated for. each., node; . In order to 
solve such a sys-tem of linear algebraic'equations, two 
basic methods exist.. The first method-is.'the-direct method, 
which is to solve the simultaneous equations'directly, 
generating an exact solution to the system of equations. 

The second me-hhod is the indirect method, which arrives at 
a solution by successive improvements of an approximate 
solution. This process is repeated until, a; desired conver¬ 
gence criteria is satisfied.. 

The basic me-thod of solution used-in thisstudy was the 
direct method.. Using- -this: method,' the-number:of calcula¬ 
tions required to achieve a solution is known, which is not 
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true af t&e iir(f±r.ec:t: rae±Jioii-. The ddrEcttiroethodi is also well 


adapted to; machiiTeE cTslciil'artions » lir. additionr,, the Gaussian 
elimin.aticnx tgchirxquES cdE sccTution , by;.'thee direet'method is 
the mast etfiicienh: methoii off salvingr aa sy;stemrof 'simultane¬ 
ous Linear aLgehraiic; ea:piat±ons. Using-this a method, however, 
extensive computer- storage, is required:.. Forrexample, a 
finite difference; aeluti-on for a gridf involving 10 nodes in 
each direction geureraEtee ro;0 simultaneouss equations for 
each iraJfndtfn... Ln. order- tcc soLve thisa probleEr by the direct 
method,, a LQCf by- LOCCC matrix: must bei stored:'for: each unknown. 

on. the basis of this V7ork, the: SDiiition- of: nonlinear 
algebraic equations using- the iterative: method proposed 
earlier is not feasible when the valuer of:the nonlinear 
terms becomes large.. For-the solutiDn-.of-J such a problem, 
the indirect method' of' solution is■ recommended. 

One possible method' of solving- such' a: system of nonlinear 
algebraic equations is ta employ a numerical! relaxation 
scheme. This method is discussed in [12]!. One advantage 
of this method is that once the relaxation pattern is 
established, the need for storing matrices is eliminated. 

Another possible method of solution is to recast the 
problem in the form of vorticity transport 


to = 
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where u is the varticitY’.. lUteE finite; diffErEBcec forms rofr 
these equations can he soilv.ext ffcor t±ie; vaiiiee ocff tfaeeunknov 7 n 
at a node in terms of the. v-seXuesE ocf: otheer unkhownssati 
surrounding nodes. FOrr instancB:,, eqnatdbnr. {22%) canr.be^: 
solved for oj in terms of at: surroundingr points; . The^ 
values obtained from the previcois iteratioir.canr.befused; 
to calculate new values.. By' chonsingr some::appropriate- 
starting pointy an iterative sodutioor. sohemee canr.bee 
established. 

A further recommendatiarr. is that: aa gradedcnetwork 'be 
used in the numerical salution.. This: graded; network should 
have finer divisions near" the jet location.. Thisrwouid: 
improve the accuracy of the. method of'modeling the-jet, 
and also provide more grid' points'near-the'jet^ , where. the 
maximum changes in flow and' temperature are. experienced. . 
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APPENDIX A 
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COMPUTER PROGRAM USED TO DEVELOP 
- FINITE DIFFERENCE OPERATORS 
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SUBROUTINE INVERT(N,EP»A,X,KERtNACT) 
IMPLICIT REAL=<--8 (A-H,0-Z) 

DIMENS ION A(NACT,NACT),X(NACT,NACT) 
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